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We study the breakpoint region on the outermost branch of current-voltage characteristics of the 
stacks with different number of intrinsic Josephson junctions. We show that at periodic boundary 
conditions the breakpoint region is absent for stacks with even number of junctions. For stacks 
with odd number of junctions and for stacks with nonperiodic boundary conditions the breakpoint 
current is increased with number of junctions and saturated at the value corresponding to the 
periodic boundary conditions. The region of saturation and the saturated value depend on the 
coupling between junctions. We explain the results by the parametric resonance at the breakpoint 
and excitation of the longitudinal plasma wave by Josephson oscillations. A way for the diagnostics 
of the junctions in the stack is proposed. 



I. INTRODUCTION 

A series of experiments devoted to intrinsic Joseph- 
son junctions 1 (UJ) shows the growing interest in 
the current-voltage characteristics (IVC) of the finite 
stacki^ii Different kinds of couplings between intrinsic 
Josephson junctions determine a variety of the IVC ob- 
served in HTSC and different models are exploited for 
their descriptions. Among them are the inductive^ and 
capacitive coupling 6,7 models, which usually give approx- 
imately similar results. A unified theory for magnetic and 
electric coupling in multistacked Josephson junctions was 
developed^ The capacitively coupled Josephson junction 
model (CCJJ) seizes the main dynamical properties of 
the IJJ system, describes the multibranch structure in 
the IVC of the stack of I J J&2 and explains the microwave 
resonant absorption^ On the other side a diffusion cur- 
rent plays an important role in the stack of IJJii. The 
CCJJ model with diffusion current (CCJJ+DC model) 
was derived in Refill on the microscopic level. It gives 
an equidistant branch structured Close to the hysteresis 
jump the system of IJJ is unstable towards a switching. 
A resonance between the Josephson and plasma oscil- 
lations causes the system to switch to another branch. 
This mechanism in case of one Josephson junction has 
been considered long time agod 

The "breakpoint region" (BPR) on the IVC of the 
stack of IJJ was demonstrated in Reffl5l and it is ex- 
plained as a result of resonance between Josephson and 
plasma oscillations. We consider, that simulation of the 
IVC of IJJ by different groups using different models 
shows the BPR on the outermost branch as well, but 
the authors did not mention it (see particularly, Fig. 3a 
in RefJjl in CCJJ model; Fig.l in Ref ill! in charge imbal- 
ance (CIB) model; Fig.2(left) in RefJjJ in CIB model). 
To our knowledge, no precise experiment to observe the 
BPR has been carried out yet. 

In this paper we show that a detailed investigation 
of the breakpoint current Ib p and BPR width Wb p gives 



us new important information concerning the creation of 
longitudinal plasma waves (LPW) in stacks of IJJ and 
the peculiarities of the stacks with a finite number of 
IJJ. We study the IVC of IJJ in the framework of the 
CCJJ+DC model and investigate the dependence of the 
Ibp and the Wb p on the number of junctions N in the 
stack. We demonstrate the existence of the BPR on the 
branches corresponding to the stack with one oscillating 
junction (O-junction) and show that such information 
may allow to develop a new method for the diagnostics 
of the IJJ. 



II. MODEL AND METHOD 

In the CCJJ model 6 - a relation between the charge 
pi and the generalized scalar potential of the Z-th 
layer is pi = — 1 5 where rrj is Debye screening 
length and is expressed through a scalar potential 
4>i and derivative of phase 9i of superconducting con- 
densate by $i{t) = cj)i - A 9®±&±L The last relation re- 
flects a nonequilibrium nature of the ac Josephson effect 
in layered superconductors 11 . The superconducting lay- 
ers are in nonequilibrium state due to the injection of 
quasiparticles and Cooper pairs. In the equilibrium state 
= 0. In the CCJJ+DC model^ with diffusion cur- 
rent J l D = — $i between layers I and I + 1 the total 
external current through the stack has a form 

j = cf + ^ + JL^ (i) 

where Vi is the voltage between superconducting layers 
/ + 1 and I (see below), tpi is the gauge-invariant phase 

difference ipi(t) = 6i+i{t) - 0i(t) - f dzA z {z, t) be- 
tween layers I + 1 and I, R is junction's resistance, A z 
is the vector potential in the barrier. This total external 
current is different from the current in the CCJJ model 
by third term in the right hand side of the equation |T|). 



2 



In the CCJJ model it is equal to Vi/R. As a result, in 
the CCJJ+DC model we obtain the following system of 
dynamical equations for the phase differences ipi 



d 2 w/dt 2 = Mv [I ~ sin ip v - Pdtpv/dt] 



(2) 



with matrix A 



/l + aG -a 
—a 1 + 2a —a 



A = 











-a 1 + 2a -a 



V 



-a 1 



(3) 



xG) 



where V runs over all N junctions, parameter a gives the 
coupling between junctions, is dissipation parameter 
(/3 2 = 1//3 C , where f3 c — u 2 R 2 C 2 is McCumber param- 
eter, Lu p is the plasma frequency and C is the capacity 
of the junction), / is external current normalized to the 
critical current I c , G = 1 + 7, 7 = s/sq = s/sjv and 
s, so, sn are the thickness of the middle, first and last 
^-layers, respectively. In the equation ((2]) time is nor- 
malized to the plasma frequency cu p r^ According to the 
proximity effect we consider that the thickness of the first 
and last layers are different from the layers inside of the 
stack. Nonperiodic boundary conditions (BC) are char- 
acterized by parameter 7 and the equations for the first 
and last layers in the system ([2]) are different from the 
equation for the middle S-layer^ii For periodic BC the 
matrix A has the form 



f 



V2a -a 

-a l + 2a -a ... 

-a 1 + 2a -a 



V 



(4) 



-a l + 2a/ 



We solve this system for the stacks with different num- 
ber N of intrinsic junctions. The numerical procedure 
has been done as follows. For a given set of model pa- 
rameters N, a, (3, 7 we simulate the IVC of the system, 
i.e. Vi(I), increasing / from zero up and then down. A 
change in the parameters N, a, /3, 7 change the branch 
structure in the IVC essentially. Their influence on the 
IVC in the CCJJ and CCJJ+DC models was discussed 
in Refsl MIMl7i . To calculate the voltages Vi(I) in each 
point of IVC (for each value of I), we simulate the dy- 
namics of the phases tpi (t) by solving system of equations 
@ using fourth order Runge-Kutta method. After simu- 
lation of the phase dynamics we calculate the dc voltages 
on each junction as 



(•5) 



where Vi is normalized to the Vo = fuv p /(2e). The aver- 
age of voltage Vi is given by 

l r Tma * 
Vi = ^ ^ / Vidt (6) 



Trnax Train JT, 



where T m i n and T max determine the interval for the 
averaging. After completing the voltage averaging for 
current / , the current I is increased or decreased by a 
small amount 51 to calculate the voltages in the next 
point of the IVC. We use a distribution of phases and 
their derivatives achieved in the previous point of the 
IVC as an initial distribution for the current point. 

Numerical stability was checked by doubling and di- 
viding in half the temporal discretization steps Dt and 
checking the influence on the IVC. Finally we can obtain 
the total dc voltage V of the stack by 



N 



v = £^ 



(7) 



1=1 



At some current / some junction (or junctions) 
switches to the nonzero voltage state and gives some 
branch of the IVC. We plot the total IVC at different 
parameters of the problem. The details concerning the 
numerical procedure are given in Refs0H3- 

To investigate the BPR in detail, we have calculated 
the IVC for different boundary conditions for the stacks 
with different number N of IJJ from N — 3 to N — 
30. For clarity we restrict the number of curves in some 
figures. 



III. IVC FOR THE STACKS WITH 
DIFFERENT NUMBER OF IJJ. 

Result of simulation of the total branch structure in 
the IVC for the stack of 10 IJJ in the CCJJ+DC model 
by the equation © is presented in the insert of Fig. QJi. 
As we can see, this IVC demonstrates the BPR on the 
outermost branch. 

The outermost branch corresponds to the state of the 
stack with all junctions in the rotating state (R-state) 17 
and it is the upper branch in the IVC. The values of 
the breakpoint current Ib p and transition current Ij (the 
jump point to the next branch in the IVC) on the out- 
ermost branch are shown by arrows in Fig. [T£l The dis- 
tance between these two values we call as the width Wb p 
of the BPR. We have found that the breakpoint current 
Ib p and BPR width Wb p depend on the parameters a and 
/3, boundary conditions and number of junctions in the 
stack. 

Let us first describe the main features of the BPR 
which follow from the results of the simulation. As we 
can see in Fig. [TJi, at 7 = both Ib p and Ij are increased 
with N, but the increase of the Ib p is monotonic. The 
IVC of the stacks with even N have larger Wb p at small 
N . The IVC in the BPR shows a chaotic behavior and 
its width is decreased with N. There is a saturation of 
N-dependence of the Ib p at large N. 

At 7 = 1 (Fig. QJ>) these features remain unchanged 
but the value of the Wb P is decreased for all stacks, espe- 
cially for N = 4 and it is equal to zero for N = 3. This 
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FIG. 1: (Color online) (a) - IVC of the outermost branch for 
the stacks with different number N of IJJ at 7 = 0; (b) - the 
same at 7 = 1. 



change in the boundary conditions leads to the relative 
changing of the BPR width in different stacks as well. 
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FIG. 2: (Color online) IVC of the outermost branch for the 
stacks with different number of junctions at periodic BC. 



the same behavior for the Ibp and BPR width Wb p for the 
stacks with odd N as in the nonperiodic case, but for the 
stacks with even N the value of the hp does not depend 
on the N and the BPR for these stacks is absent. 

In Fig. [3^, the Ib p as a function of N for 7 = (squares, 
curve 1) and periodic (circles, curve 2) boundary condi- 
tions at different value of the coupling parameter a is pre- 
sented. We stress the coincidence of the N-dependencies 
of the Ibp for stacks with odd numbers of IJJ for peri- 
odic and 7 = cases. The increase in a leads to the 
saturation of the N-dependence at larger N. The value of 
saturated Ibp is decreased with coupling and consists of 
0.576 at a = 1, 0.454 at a = 0.5 and 0.304 at a = 0.1. At 
a = the breakpoint coincides with the return current, 
so the Ibp has the same value for the stacks with different 
number of junctions. 

The N-dependence of the BPR width Wb p for stacks 
with even and odd number of junctions at 7 = (curves 
leven and l dd) and periodic (curve 2 a dd) boundary con- 
ditions is shown in Fig. [3)d. The main feature here is a 
decrease of the BPR width Wb p with N at large N. At 
small N in the interval (3,6) we observe the increase of 
the Wbp with N. 



IV. THE ORIGIN OF THE BREAKPOINT ON 
THE OUTERMOST BRANCH 

To explain the observed features of the finite stacks 
IVC let us discuss the origin of the breakpoint on the 
outermost branch. The hysteresis jump in the IVC is 
associated with the change of the distribution pattern of 
rotating phase motions ji£ But the question "why does a 
change in the current leads to the change in a distribution 
pattern of the rotating phase motions" is still open. We 
consider a case that all junctions are in the rotating state, 



the time average of tpi (ipi 



■ \ T ma " fidt) 



1.6. U-L-LXXVj CXV^XC^g,^ Ul rp J rp 

is constant and that of sin ipi is zero for these junctions. 
For the oscillating junctions the situation is opposite: the 
time average of (pi is zero and that of sintp/ is constant. 

As we mentioned above, the outermost branch in the 
IVC corresponds to the state of the stack with all junc- 
tions in the rotating state. Let us write an equation for 
the difference of phase differences Si = tpi+i — fi for the 
outermost branch. 

By subtracting equation ((2|) for (/ + l)th and (l)th 
junctions we get 



l+i 



<Pl) 



(1 



*V (2) ){sin( W+1 ) 



- sin(ipz) 



(8) 



Here V^/z = fi+i + fi-i — 2fi is the discrete Laplacian. 
Consider a linear approximation sin(y>;+i) — sin(</?;) « 
6icos(tp), where <p ~ Qt = -hVt, is Josephson fre- 
quency, V is total voltage of the stack, we obtain 



The IVC at periodic boundary condition (Fig. ^ show 



^ + (l-aV (2) )(cos(^/+^)=0 (9) 
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FIG. 3: (Color online) (a) - the N-dependence of the h p for 
7 = (curve 1) and periodic (curve 2) boundary conditions 
at different a; (b) - the N-dependence of the BPR width Wb p 
for stacks with even and odd number of junctions at 7 = 
(1) and periodic (2) boundary conditions. 



Expanding 5i(t) in the Fourier series 

<yj(i) = 5D* fc( 



, //>•/ 



(10) 



the linearized equation for the Fourier component of a 
difference of the phase differences 5k between neighbor- 
junctions can be written in the form 15 



6k + P{k)5 k + cos(n(k)T)6 k = 0, 



(11) 



where r = u> p (k)t, Lo p (k) = u> p C a , (3(k) — f3C a , f2(fc) = 
n/C a and C a = ^l + 2a(l-cos(fc)). 

The important fact for us is that this linearized equa- 
tion manifests a parametric resonance in the system of 
IJJ. In Fig. UK we plot the resonance region for this equa- 
tion on the diagram (3(k) — fi(fe). The dark stripe on this 
figure is actually the distribution of the dots, correspond- 
ing to the positions of the breakpoints of the outermost 
branch. Using the breakpoint values of the voltage in the 
equation fi(fc) = f2/C Q = V/NC a we obtain this distri- 
bution of the breakpoints by the variation of the coupling 
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FIG. 4: (Color online) (a) - Parametric resonance region in 
Q,(k) — f3(k) diagram. Black dots (stripe) correspond to the 
breakpoint current It, p in the IVC for k = n at different values 
of parameters a and /3; (b) - Charge on the S-layers in the 
stacks with 10 IJJ (bl) and 11 IJJ (b2) at periodic BC and 
with 10 IJJ at 7 = 1 (b3). 



parameter a in the interval (1,2) with a step 0.1 and the 
dissipation parameter f3 in the interval (0.01,0.34) with 
the step 0.01 at each value of a. In contrast to the re- 
sults presented in Ref[l5| where the positions of the dots 
on the diagram (3{k) — f2(fc) were obtained by crude esti- 
mation, here we have done the precise numerical calcula- 
tions. This calculations show more close positions of the 
breakpoints to the boundary of the resonance region in 
the chosen intervals of a and (3. The reason, that the po- 
sition of the breakpoints are shifted from the boundary 
of the resonance region in Fig. [4^,, is the linear approxi- 
mation used to obtain the equation @. 

The breakpoints are inside of the resonance region, i.e. 
the resonance between the Josephson and plasma oscil- 
lations is approached at the breakpoint Ib p . As a result 
the plasma mode is excited by the Josephson oscillations. 
We can prove this statement directly. By Maxwell equa- 
tion, div(E/d) — 4-7rp, we express the charge pi on the 
S-layer I by the voltages Vij-i and Vij+i in the neighbor 
insulating layers pi = 3^3(^+1 - Vj_i,j). The time de- 
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pendence of the pi, presented in Fig. |4}d1, demonstrates 
that at periodic BC in the stacks with 10 IJJ the LPW 
with k = 7r is realized. Really, the charge on the near- 
est neighbor layers has the same value and an opposite 
sign. Fig. [4})2 shows the distribution of the charge on 
the S-layers in the stacks with 11 IJJ at periodic BC. 
In this case we observe the creation of the LPW with 
k = 107r/l 1. To determine the mode of the LPW at 
nonperiodic BC need more detail investigation. Fig. (4ja3 
demonstrates the charge distribution in the stack with 10 
IJJ at 7 = 1 near the breakpoint. 

The wavelengths of the standing LPW which can be 
realized in the stack with N junctions are N/n units of 
lattice in z-direction, where n changes from 1 to N/2 
for stacks with even number of junctions and from 1 
to (N - l)/2 for odd N. The voltage of the stack at 
the breakpoint is related to the wave number k of the 
LPW by formulae V = NQ(k) yjl + 2a(l - cos(fc)), so 
the largest breakpoint voltage V in the current decreas- 
ing process corresponds to the creation of the LPW with 
k equal to it (7r-mode) for stacks with even number of IJJ 
and modes with k = (N — 1)tt/N for stacks with odd N. 



V. DISCUSSION OF THE MAIN RESULTS 

Let us now return to the results presented in Fig. [TJ- 
Fig. [3] and demonstrate that they are in agreement with 
the ideas stated in the previous section. 

According to these ideas, in the stack with 10 IJJ the 
LPW with k — (N — 1)tt/N is created and it leads to the 
increase of the Ib p with N and its saturation to the value, 
corresponding to the 7r-mode. The same modes are cre- 
ated in the stacks with 7=1 which outermost branches 
of the IVC are presented in Fig. QJ>. They demonstrate 
the saturation of the Ibp to the same value. 

As Fig. [2] shows, at periodic BC we observe the same 
value of Ibp in the all stacks with even N. It is in agree- 
ment with our suggestion that in this case the LPW with 
k = 7r is created. We check it directly as well, by the time 
dependence of the pi. We find that at periodic BC in the 
stacks with even N the charge on the nearest neighbor 
layers has the same value and an opposite sign which 
means that the LPW with k — ir is realized. 

In the stacks with odd N the 7r-mode cannot exist, so 
as we mentioned in the previous section, the LPW with 
the largest k equal to k = (N — l)w/N is created. The 
creation of different modes of the LPW leads to the dif- 
ferent Ibp and this fact explains the difference in the IVC 
at periodic BC of the stacks with even and odd numbers 
of IJJ. With increase in N the wave number k limits to 
7r and it leads to the increase in Ibp which we observe in 
Fig.H 

It explains as well the saturation of the Ibp to the value 
corresponding to the Ibp for stacks with even N which is 
demonstrated in Fig. [3^. 

The difference between the charge distribution on the 
S-layers at the breakpoint current at fixed time for stacks 



with even and odd numbers of IJJ (10 and 11) at periodic 
BC is demonstrated in Fig. \^p. In the stack with odd N 
the charge on the first and last layers oscillates in-phase 
and the oscillations on the neighbor S-layers are different 
from the oscillations in 7r-mode. Because at 7 = the 
charge on the first and last layers is screened due to the 
proximity effect, the LPW with the k = (N — \)ir/N 
as at periodic BC for stacks with odd N is also created 
both for stacks with odd and even N. This is a reason 
the values of Ib p for stacks with odd N coincide for both 
these BC as shown in Fig. 

At 7 = 1 the same modes are created, but the character 
of the charge distribution among the S-layers is different 
(as we can see in Fig. \^p3 ) . As a result in this case the 
N-dependence of the Ibp is stronger than at 7 = 0, but at 
N — > 00 it is saturated as well and the saturation value 
is the Ibp for even junctions at periodic BC. 

The influence of the coupling parameter on the value of 
the breakpoint current Ib p and its N-dependence which 
is demonstrated in Fig. [3^ has a clear explanation. It 
consists in the mentioned in the previous section the 
a-dependence of the voltage, and correspondingly, the 
breakpoint current Ib p - They are proportional to the 
term yl + 2a(l — cos(fc)). An increase in k leads to the 
term proportional to the yT+4a. The decrease in a 
makes this influence weaker. 

We may explain an increase in Wb P at small N, which 
is shown in Fig. [HJd by a commensurability effect on the 
width of the BPR. At periodic BC Wb p = 0, if the wave- 
length A of the LPW is A — n, where n — 2,3, 4... lattice 
units in z-direction. As we can see in Fig. [21 at A = 3 (for 
stacks with N = 3) and A = 2 (7r-mode), the IVC do not 
show the BPR. The creation of the LPW with wavelength 
in the middle of the interval (2,3) should correspond to 
some maximum of the BPR width Wb P . In the stack with 
5 IJJ the LPW with A = 2.5 (k = 4tt/5) is created. It 
does not coincide exactly with the result obtained by the 
simulation, because we use for the explanation the lin- 
earized equation for difference of phase differences and 
Fourier expansion in the finite stacks. The wave number 
k is not well defined in this case. With increase in N the 
wave number k of LPW limits to it and it explains the 
increase in Ibp and the decrease in Wb p at large N which 
demonstrated by Fig. [3}x The change in the boundary 
conditions changes the character of the charge oscilla- 
tions on the S-layers. Particularly, at 7 = 1 we observe 
a decrease in the Wb P for all stacks. Our analysis shows 
that in stacks with even N at nonperiodic BC the charge 
on the second and N — 1 layers oscillate in-phase, but 
anti-phase in the stacks with odd N. We consider that 
the value of the BPR width Wb p depends on the charac- 
ter of the charge oscillations on the S-layers. 

Finally, we note that in the case of coupling between 
junctions the parameter [3 cannot be determined in a 
usual way by the return current, because it depends now 
on two parameters, and a. The dependence of the 
Ibp and BPR width Wb p on the dissipation and coupling 
parameters opens an opportunity to develop the new 
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method for determination of these parameters for the 
stacks of IJJ. This question will be discussed in detail 
somewhere else. 



VI. ONE OSCILLATING JUNCTION 

Let us now discuss briefly the breakpoints on the other 
branches of the IVC. 
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FIG. 5: (Color online) The BPR on the branches of the IVC 
of the stacks with one oscillating junction in the case of 10 IJJ 
at a = 1, P = 0.2 and 7 = 0. The upper curve corresponds 
to the real scale of voltage, but the other ones are shifted for 
clarity by two units down. 

As we mentioned above a resistive state in the sys- 
tem of IJJ is realized as a state with different number of 
R- and O-junctions.?*^ The different positions of R- and 
O-junctions in the stack (different configurations) corre- 
spond to the different states of IJJ system. The equidis- 
tant positions of the O-junction from the ends of the stack 
(for example, the stacks with 1st or 10th O-junction) lead 
to the same state. So, there are five different states in the 
stack with one O-junction corresponding to the different 
position of this junction. Fig. \5\ shows the BPR on the 
branches of the IVC of the stacks with one O-junction in 
the case of 10 IJJ at a = 1, = 0.2 and 7 = 0. 

The equidistant positions of the O-junction from the 
ends of the stack lead to the same value of Ib p and the 
same width of BPR. The shift of the O-junction from 
the end of the stack to its center decreases the 7& p of the 
corresponding state. So, we may establish a delay of the 
LPW creation in the current decreasing process when the 
position of the O-junction is shifted to the center of the 
stack. 

We consider that the origin of such behaviour is the 
next. This one oscillating junction separates the stack 
in two parts with different number of R-j unctions which 
are weakly coupled through it. With a decrease in cur- 



rent the first LPW is created in the part with the largest 
Ibp (with the largest number of junctions). The shift 
of the O-junction and the decrease in the number of R- 
junctions in this part lead to the decrease of the Ib p as 
Fig. [T] demonstrates. The increase of the number of junc- 
tions in the second part might manifest the second break- 
point which is related to the creation of LPW in this 
second part of the stack. Such situation is observed for 
N = 10 when O-junction occupies the 5th or 6th site in 
the stack. The width of the BPR in the other branches 
of IVC depends essentially on the state of the stack. 

For the other branches the increase in the number of 
the O-junctions in the stack decreases the number of ef- 
fective junctions for creation of the LPW and it leads 
to the decrease of the return current. This fact explains 
why we can obtain a total branch structure in the hys- 
teresis region, because in the other case we would not be 
able to observe it in the simulation. The correspondence 
between the position of the O-junction in the stack and 
the value of the Ib p opens the possibility for junction di- 
agnostics, i.e. by measuring the value of the Ib p we can 
answer the question which junction in the stack goes into 
R- or O-state. From the other side, the monitoring of the 
transitions between branches is useful for understanding 
the phase dynamics in the system of IJJ. 



VII. CONCLUSIONS 

In conclusion, we stress that the BPR in the IVC "nat- 
urally" follows from the solution of the system of the dy- 
namical equations for the phase difference for the stack of 
IJJ. In breakpoint region the plasma mode is a station- 
ary solution of the system and this fact might be used in 
some applications, particularly, in high frequency devices 
such as THz oscillators and mixers. The detailed study of 
the breakpoint current and breakpoint region width gives 
a new opportunity for the investigation of the properties 
of IJJ and and develop new methods for the determina- 
tion of the parameters of IJJ and diagnostic of IJJ in the 
stacks. 
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